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Abstract 


We describe an embarrassingly parallel, anytime Monte Carlo method for 
likelihood-free models. The algorithm starts with the view that the stochastic- 
ity of the pseudo-samples generated by the simulator can be controlled externally 
by a vector of random numbers u , in such a way that the outcome, knowing u , 
is deterministic. For each instantiation of u we run an optimization procedure to 
minimize the distance between summary statistics of the simulator and the data. 
After reweighing these samples using the prior and the Jacobian (accounting for 
the change of volume in transforming from the space of summary statistics to the 
space of parameters) we show that this weighted ensemble represents a Monte 
Carlo estimate of the posterior distribution. The procedure can be run embar¬ 
rassingly parallel (each node handling one sample) and anytime (by allocating 
resources to the worst performing sample). The procedure is validated on six ex¬ 
periments. 


1 Introduction 

Computationally demanding simulators are used across the full spectrum of scientific and industrial 
applications, whether one studies embryonic morphogenesis in biology, tumor growth in cancer 
research, colliding galaxies in astronomy, weather forecasting in meteorology, climate changes in 
the environmental science, earthquakes in seismology, market movement in economics, turbulence 
in physics, brain functioning in neuroscience, or fabrication processes in industry. Approximate 
Bayesian computation (ABC) forms a large class algorithms that aims to sample from the posterior 
distribution over parameters for these likelihood-free (a.k.a. simulator based) models. Likelihood- 
free inference, however, is notoriously inefficient in terms of the number of simulation calls per 
independent sample. Further, like regular Bayesian inference algorithms, care must be taken so that 
posterior sampling targets the correct distribution. 

The simplest ABC algorithm, ABC rejection sampling, can be fully parallelized by running indepen¬ 
dent processes with no communication or synchronization requirements. I.e. it is an embarrassingly 
parallel algorithm. Unfortunately, as the most inefficient ABC algorithm, the benefits of this ti¬ 
tle are limited. There has been considerable progress in distributed MCMC algorithms aimed at 
large-scale data problems GJH1. Recently, a sequential Monte Carlo (SMC) algorithm called “the 
particle cascade” was introduced that emits streams of samples asynchronously with minimal mem¬ 
ory management and communication ini. In this paper we present an alternative embarrassingly 
parallel sampling approach: each processor works independently, at full capacity, and will indefi¬ 
nitely emit independent samples. The main trick is to pull random number generation outside of the 
simulator and treat the simulator as a deterministic piece of code. We then minimize the difference 
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between observations and the simulator output over its input parameters and weight the final (opti¬ 
mized) parameter value with the prior and the (inverse of the) Jacobian. We show that the resulting 
weighted ensemble represents a Monte Carlo estimate of the posterior. Moreover, we argue that the 
error of this procedure is 0(e) if the optimization gets e-close to the optimal value. This “Opti¬ 
mization Monte Carlo” (OMC) has several advantages: 1) it can be run embarrassingly parallel, 2) 
the procedure generates independent samples and 3) the core procedure is now optimization rather 
than MCMC. Indeed, optimization as part of a likelihood-free inference procedure has recently been 
proposed IH; using a probabilistic model of the mapping from parameters to differences between 
observations and simulator outputs, they apply “Bayesian Optimization” (e.g. USED) to efficiently 
perform posterior inference. Note also that since random numbers have been separated out from the 
simulator, powerful tools such as “automatic differentiation” (e.g. (H) are within reach to assist 
with the optimization. In practice we find that OMC uses far fewer simulations per sample than 
alternative ABC algorithms. 

The approach of controlling randomness as part of an inference procedure is also found in a related 
class of parameter estimation algorithms called indirect inference fTTTl . Connections between ABC 
and indirect inference have been made previously by Q as a novel way of creating summary statis¬ 
tics. An indirect inference perspective led to an independently developed version of OMC called the 
“reverse sampler” mm. 

In Section [2] we briefly introduce ABC and present it from a novel viewpoint in terms of random 
numbers. In Section [5] we derive ABC through optimization from a geometric point of view, then 
proceed to generalize it to higher dimensions. We show in Section [4] extensive evidence of the 
correctness and efficiency of our approach. In Section [5] we describe the outlook for optimization- 
based ABC. 

2 ABC Sampling Algorithms 

The primary interest in ABC is the posterior of simulator parameters 6 given a vector of (statistics 
of) observations y, p(0 |y). The likelihood p(y\0) is generally not available in ABC. Instead we can 
use the simulator as a generator of pseudo-samples x that reside in the same space as y. By treating 
x as auxiliary variables, we can continue with the Bayesian treatment: 

/a i n p(0)p(y|0) _ p(fl)/p e (y|x)p(x|6>)dx 

p( y) ~ /p(0)/pe(y|x)p(x|0) dxdd 

Of particular importance is the choice of kernel measuring the discrepancy between observations y 
and pseudo-data x. Popular choices for kernels are the Gaussian kernel and the uniform e-tube/ball. 
The bandwidth parameter e (which may be a vector e accounting for relative importance of each 
statistic) plays critical role: small e produces more accurate posteriors, but is more computationally 
demanding, whereas large e induces larger error but is cheaper. 

We focus our attention on population-based ABC samplers, which include rejection sampling, im¬ 
portance sampling (IS), sequential Monte Carlo (SMC) I6ll20l and population Monte Carlo 0. In 
rejection sampling, we draw parameters from the prior 0 ~ p(0), then run a simulation at those 
parameters x ~ p(x|0); if the discrepancy p(x, y) < e, then the particle is accepted, otherwise it 
is rejected. This is repeated until n particles are accepted. Importance sampling generalizes rejec¬ 
tion sampling using a proposal distribution q^O) instead of the prior, and produces samples with 
weights Wi oc p(0)/q{0). SMC extends IS to multiple rounds with decreasing e, adapting their par¬ 
ticles after each round, such that each new population improves the approximation to the posterior. 
Our algorithm has similar qualities to SMC since we generate a population of n weighted particles, 
but differs significantly since our particles are produced by independent optimization procedures, 
making it completely parallel. 

3 A Parallel and Efficient ABC Sampling Algorithm 

Inherent in our assumptions about the simulator is that internally there are calls to a random number 
generator which produces the stochasticity of the pseudo-samples. We will assume for the moment 
that this can be represented by a vector of uniform random numbers u which, if known, would 
make the simulator deterministic. More concretely, we assume that any simulation output x can 
be represented as a deterministic function of parameters 0 and a vector of random numbers u , 
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(a) Do = D y 


(b) D e < D y 


Figure 1: Illustration of OMC geometry, (a) Dashed lines indicate contours f(0 : u) over 6 for several u. For 
three values of it, their initial and optimal 6 positions are shown (solid blue/white circles). Within the grey 
acceptance region, the Jacobian, indicated by the blue diagonal line, describes the relative change in volume 
induced in f(0 : it) from a small change in 9. Corresponding weights oc 1/| J\ are shown as vertical stems, (b) 
When Do < D y , here 1 < 2, the change in volume is proportional to the length of the line segment inside the 
ellipsoid (| J T J| 1 ' 2 ). The orange line indicates the projection of the observation onto the contour of f(0 , it) (in 
this case, identical to the optimal). 


i.e. x = This assumption has been used previously in ABC, first in “coupled ABC” 

El and also in an application of Hamiltonian dynamics to ABC m. We do not make any further 
assumptions regarding u or p(u), though for some problems their dimension and distribution may be 
known a priori. In these cases it may be worth employing Sobol or other low-discrepancy sequences 
to further improve the accuracy of any Monte Carlo estimates. 

We will first derive a dual representation for the ABC likelihood function p e (y\0) (see also ITbl ). 

Pe(y|0) = J Pe(y|x)p(x|0) dx = y , y'p e (y|x)I[x = f(0,u)]p(u) dxdu (2) 

= jPe(y\f(0,u))p(u) du (3) 

leading to the following Monte Carlo approximation of the ABC posterior, 

Pe(0 ly) oc p(0) Jp(u)p e (y\f(u, 6)) du sa J ^p e (y|/ («,, 0))p{G) u t ~ p(u) (4) 

i 

Since p e is a kernel that only accepts arguments y and f(ui, 6) that are e close to each other (for 
values of e that are as small as possible), Equation [4] tells us that we should first sample values for 
u from p(u) and then for each such sample find the value for 0° that results in y = f(0°, u). In 
practice we want to drive these values as close to each other as possible through optimization and 
accept an 0(e) error if the remaining distance is still 0(e). Note that apart from sampling the values 
for u this procedure is deterministic and can be executed completely in parallel, i.e. without any 
communication. In the following we will assume a single observation vector y, but the approach is 
equally applicable to a dataset of N cases. 

3.1 The case Dq = D y 

We will first study the case when the number of parameters 6 is equal to the number of summary 
statistics y. To understand the derivation it helps to look at Figure [la| which illustrates the derivation 
for the one dimensional case. In the following we use the following abbreviation: fi(0) stands for 
f(0,Ui). The general idea is that we want to write the approximation to the posterior as a mixture 
of small uniform balls (or delta peaks in the limit): 

p(o |y) ~ ~y2pe(y\f(ui,e))p(0 )« - 'y^w i u € (o\e*)p(o) (5) 
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with Wi some weights that we will derive shortly. Then, if we make e small enough we can replace 
any average of a sufficiently smooth function h{0) w.r.t. this approximate posterior simply by eval¬ 
uating h(0) at some arbitrarily chosen points inside these balls (for instance we can take the center 
of the ball 0*), 

[ h(0)p(6 |y) dd « -^h(d*) w t p(0*) ( 6 ) 

J n i 

To derive this expression we first assume that: 

Pe(y|/ i W) = C(e)I[||y-/ i (0)|| 2 <e 2 ] (7) 

i.e. a ball of radius e. C(e) is the normalizer which is immaterial because it cancels in the posterior. 
For small enough e we claim that we can linearize fi(0) around 0°: 

f i {e) = f i {e°)+r i {e-o°) + Ri, r { = o(\\o - e?\\ 2 ) ( 8 ) 

where J° is the Jacobian matrix with columns - . We take 6° to be the end result of our 
optimization procedure for sample Ui. Using this we thus get, 

lly - fi(m 2 ~ ll(y - fiW)) - j tie - o?) - iu || 2 (9) 

We first note that since we assume that our optimization has ended up somewhere inside the ball 
defined by ||y — fi(0 )|| 2 < e 2 we can assume that ||y — fi(0°)\\ = 0(e). Also, since we only 
consider values for 0 that satisfy ||y — fi(0)\\ 2 < e 2 , and furthermore assume that the function 
fi(0) is Lipschitz continuous in 0 it follows that ||0 — 0°\\ = 0(e) as well. All of this implies 
that we can safely ignore the remaining term Ri (which is of order 0(\ \0 — 0°\ | 2 ) = 0(e 2 )) if we 
restrict ourselves to the volume inside the ball. 

The next step is to view the term I[||y — fi(0)\\ 2 < e 2 ] as a distribution in 0. With the Taylor 
expansion this results in, 

i[(0 - % - Ji ,_1 (y - /i(fl?))) T Ji r J?(0 - % - J°’ _1 (y - /<(*?))) < e 2 ] (io) 

This represents an ellipse in 0-space with a centroid 0* and volume Vi given by 

e* = d° + i°r\y - MW) ^ = 1 (ii) 

y det(J° T J°) 

with 7 a constant independent of i. We can approximate the posterior now as, 

K i \/det(jfj“) K i dtetOfK) 

where in the last step we have send e 0. Finally, we can compute the constant k through nor¬ 
malization, k, = £.p(0*) det(J^ T J^) -1 / 2 . The whole procedure is accurate up to errors of the 
order 0(e 2 ), and it is assumed that the optimization procedure delivers a solution that is located 
within the epsilon ball. If one of the optimizations for a certain sample U{ did not end up within 
the epsilon ball there can be two reasons: 1 ) the optimization did not converge to the optimal value 
for 0 , or 2 ) for this value of u there is no solution for which f(0\u) can get within a distance e 
from the observation y. If we interpret e as our uncertainty in the observation y, and we assume that 
our optimization succeeded in finding the best possible value for 0 , then we should simply reject 
this sample 0^. However, it is hard to detect if our optimization succeeded and we may therefore 
sometimes reject samples that should not have been rejected. Thus, one should be careful not to 
create a bias against samples Ui for which the optimization is difficult. This situation is similar to a 
sampler that will not mix to remote local optima in the posterior distribution. 

3.2 The case Dq < D y 

This is the overdetermined case and here the situation as depicted in Figure [lb] is typical: the mani¬ 
fold that f(0,Ui) traces out as we vary 0 forms a lower dimensional surface in the D y dimensional 
enveloping space. This manifold may or may not intersect with the sphere centered at the observa¬ 
tion y (or ellipsoid, for the general case e instead of e). Assume that the manifold does intersect the 
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epsilon ball but not y. Since we trust our observation up to distance e, we may simple choose to 
pick the closest point 0 * to y on the manifold, which is given by, 

0* = 0° + jf(y - jf = (13) 

where J°^ is the pseudo-inverse. We can now define our ellipse around this point, shifting the 
center of the ball from y to (which do not coincide in this case). The uniform distribution 

on the ellipse in 0-space is now defined in the Dq dimensional manifold and has volume Vi = 
7 det(J^ T J^) -1 / 2 . So once again we arrive at almost the same equation as before (Eq.|l2|) but with 
the slightly different definition of the point 0* given by Eq.[l3] Crucially, since | |y — m ) II <e 2 
and if we assume that our optimization succeeded, we will only make mistakes of order 0{e 2 ). 

3.3 The case Dq > D y 

This is the underdetermined case in which it is typical that entire manifolds (e.g. hyperplanes) may 
be a solution to ||y — fi(0*)\\ = 0. In this case we can not approximate the posterior with a 
mixture of point masses and thus the procedure does not apply. However, the case Dq > D y is less 
interesting than the other ones above as we expect to have more summary statistics than parameters 
for most problems. 

4 Experiments 

The goal of these experiments is to demonstrate 1) the correctness of OMC and 2) the relative 
efficiency of OMC in relation to two sequential MC algorithms, SMC (aka population MC ED and 
adaptive weighted SMC s To demonstrate correctness, we show histograms of weighted samples 
along with the true posterior (when known) and, for three experiments, the exact OMC weighted 
samples (when the exact Jacobian and optimal 0 is known). To demonstrate efficiency, we compute 
the mean simulations per sample (SS)—the number of simulations required to reach an e threshold— 
and the effective sample size (ESS), defined as l/w T w. Additionally, we may measure ESS/n, the 
fraction of effective samples in the population. ESS is a good way of detecting whether the posterior 
is dominated by a few particles and/or how many particles achieve discrepancy less than epsilon. 

There are several algorithmic options for OMC. The most obvious is to spawn independent pro¬ 
cesses, draw u for each, and optimize until e is reached (or a max nbr of simulations run), then 
compute Jacobians and particle weights. Variations could include keeping a sorted list of discrepan¬ 
cies and allocating computational resources to the worst particle. However, to compare OMC with 
SMC, in this paper we use a sequential version of OMC that mimics the epsilon rounds of SMC. 
Each simulator uses different optimization procedures, including Newton’s method for smooth sim¬ 
ulators, and random walk optimization for others; Jacobians were computed using one-sided finite 
differences. To limit computational expense we placed a max of 1000 simulations per sample per 
round for all algorithms. Unless otherwise noted, we used n = 5000 and repeated runs 5 times; lack 
of error bars indicate very low deviations across runs. We also break some of the notational conven¬ 
tion used thus far so that we can specify exactly how the random numbers translate into pseudo-data 
and the pseudo-data into statistics. This is clarified for each example. Results are explained in 
Figures [4] to [4] 

4.1 Normal with Unknown Mean and Known Variance 

The simplest example is the inference of the mean 0 of a univariate normal distribution with known 
variance a 2 . The prior distribution 7 r(0) is normal with mean $o and variance ka 2 , where k > 0 is 
a factor relating the dispersions of 0 and the data y n . The simulator can generate data according to 
the normal distribution, or deterministically if the random effects r Urn are known: 

K(Xm\0) = N(Xm\0',(J 2 ) => x rn =e + r Urn (14) 

where r Um = ay/2 erf _1 (2u rn — 1) (using the inverse CDF). A sufficient statistic for this problem is 
the average s(x) = ^ x Therefore we have f(0, u) = 0 + R{u) where R(u) = r um/M 
(the average of the random effects). In our experiment we set M = 2 and y = 0. The exact 
Jacobian and 6° can be computed for this problem: for a draw m, Ji = 1; if s( y) is the mean of the 
observations y, then by setting f{0°,Ui) = s( y) we find 0° = s( y) — R(u{). Therefore the exact 
weights are Wi oc 7r( 9 ?), allowing us to compare directly with an exact posterior based on our dual 
representation by u (shown by orange circles in Figure [2] top-left). We used Newton’s method to 
optimize each particle. 
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Unknown Mean (eps = 0.01) Mixture of Normals (eps = 0.01) 





Figure 2: Left: Inference of unknown mean. For e 0.1, OMC uses 3.7 SS; AW/SMC uses 20/20 SS; at e 0.01, 
OMC uses 4 SS (only 0.3 SS more), and SMC jumps to 110 SS. For all algorithms and e values, ESS/n=l. 
Right: Inference for mixture of normals. Similar results for OMC; at e 0.025 AW/SMC had 40/50 SS and at 
e 0.01 has 105/120 SS. The ESS/n remained at 1 for OMC, but decreased to 0.06/0.47 (AW/SMC) at e 0.025, 
and 0.35 for both at e 0.01. Not only does the ESS remain high for OMC, but it also represents the tails of the 
distribution well, even at low e. 


4.2 Normal Mixture 

A standard illustrative ABC problem is the inference of the mean 0 of a mixture of two normals 
|[T9l l3ll5ll: p(x\0) = pN(0, (j\) + (1 — p) J\f(0, a 2 ), with tt(0) = U(6 a , 0b) where hyperparameters 
are p = 1/2, a 2 = 1, a 2 = 1/100, 0 a = —10, 6b = 10, and a single observation scalar y = 0. 

For this problem M = 1 so we drop the subscript m. The true posterior is simply p(6\y = 0) oc 

p A f{0, cr 2 ) + (1 — p) A/*(#, erf), 0 G { — 10,10}. In this problem there are two random numbers u\ 
and u 2 , one for selecting the mixture component and the other for the random innovation; further, 
the statistic is the identity, i.e. s(x ) = x: 

x = [ui < p\(0 + <Ji\/2erf(2tx 2 - 1)) + [ui > p\(0 + cr 2 \/2erf(2tx 2 - 1)) (15) 

= 0 + a/ 2 erf(2tx 2 — l)cr[ Ul< ^cr 2 Ul “^ = 6 + R(u) (16) 

where R(u) = a/ 2 erf (2 u 2 — l)cr[ Ul< ^ cr^ 1 -^. As with the previous example, the Jacobian is 1 
and 0° = y — R{ui) is known exactly. This problem is notable for causing performance issues in 
ABC-MCMC Q3 and its difficulty in targeting the tails of the posterior 0; this is not the case for 
OMC. 


4.3 Exponential with Unknown Rate 

In this example, the goal is to infer the rate 0 of an exponential distribution, with a gamma prior 
p(0) = Gamma(#|a, /3), based on M draws from Exp($): 

p(x m \0) = Exp(x m |6>) = 6 ex.p(-6xm) x m = - ^ ln(l - u m ) = ^ r Um (17) 

where r Um = — ln(l — u m ) (the inverse CDF of the exponential). A sufficient statistic for this 
problem is the average s(x) = . Again, we have exact expressions for the Jacobian 

and 0°, using f(0,Ui) = R(ui)/0, Ji = —R(ui)/6 2 and 6° = R(ui)/s(y). We used M = 2, 
s( y) = 10 in our experiments. 

4.4 Linked Mean and Variance of Normal 

In this example we link together the mean and variance of the data generating function as follows: 

p(x m \6) = Af(x m \0,6 2 ) => x m = 6 + 0V2 erf -1 (2u m - 1) = 0r Um (18) 
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Figure 3: Left: Inference of rate of exponential. A similar result wrt SS occurs for this experiment: at e 1, 
OMC had 15 v 45/50 for AW/SMC; at e 0.01, SS was 28 OMC v 220 AW/SMC. ESS/n dropping with below 
1: OMC drops at e 1 to 0.71 v 0.97 for SMC; at e 0.1 ESS/n remains the same. Right: Inference of linked 
normal. ESS/n drops significantly for OMC: at e 0.25 to 0.32 and at e 0.1 to 0.13, while it remains high 
for SMC (0.91 to 0.83). This is the result the inability of every Ui to achieve p < e, whereas for SMC, the 
algorithm allows them to “drop” their random numbers and effectively switch to another. This was verified 
by running an expensive fine-grained optimization, resulting in 32.9% and 13.6% optimized particles having 
p under e 0.25/0.1. Taking this inefficiency into account, OMC still requires 130 simulations per effective 
sample v 165 for SMC (ie 17/0.13 and 136/0.83). 


where r Um = 1 + y/2 erf 1 (2u m — 1). We put a positive constraint on 0: p(6) = U( 0,10). We used 
2 statistics, the mean and variance of M draws from the simulator: 

Sl(x) = JjXm => fi{0,u) = 0R(u) d ^QQ — = R ( u ) ( 19 ) 

S 2 ( x ) = T - s a ( x )) 2 =>- f 2 (0,u) = e 2 V(u) — = 26V(u) (20) 

m 

where V(u ) = ^Z m ^ rn /M — R(u ) 2 and R(u ) = '}2 m r Urn /M\ the exact Jacobian is therefore 
[. R(u ), 2 0V(u)] T . In our experiments M = 10, s( y) = [2.7,12.8]. 


4.5 Lotka-Volterra 

The simplest Lotka-Volterra model explains predator-prey populations over time, controlled by a set 
of stochastic differential equations: 

^ = 6 iXi - e 2 x 1X2 + T\ ^ = -0 2 x 2 - 0 3 x !X 2 + r 2 (21) 

at at 

where x\ and x 2 are the prey and predator population sizes, respectively. Gaussian noise r ~ 
A/"(0,10 2 ) is added at each full time-step. Lognormal priors are placed over 6. The simulator 
runs for T = 50 time steps, with constant initial populations of 100 for both prey and predator. 
There is therefore P = 2T outputs (prey and predator populations concatenated), which we use 
as the statistics. To run a deterministic simulation, we draw Ui ~ 7r(u) where the dimension of 
u is P. Half of the random variables are used for r\ and the other half for r 2 . In other words, 
r Ust = IOa/ 2 erf _1 (2ti s t — 1), where s G {1, 2} for the appropriate population. The Jacobian is a 
100 x 3 matrix that can be computed using one-sided finite-differences using 3 forward simulations. 


4.6 Bayesian Inference of the M/G/l Queue Model 

Bayesian inference of the M/G/l queuing model is challenging, requiring ABC algorithms 00 or 
sophisticated MCMC-based procedures m. Though simple to simulate, the output can be quite 
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Figure 4: Top: Lotka-Volterra. Bottom: M/G/l Queue. The left plots shows the posterior mean ±1 std errors 
of the posterior predictive distribution (sorted for M/G/l). Simulations per sample and the posterior of 0\ are 
shown in the other plots. For L-V, at e 3, the SS for OMC were 15 v 116/159 for AW/SMC, and increased 
at e 2 to 23 v 279/371. However, the ESS/n was lower for OMC: at e 3 it was 0.25 and down to 0.1 at e 2, 
whereas ESS/n stayed around 0.9 for AW/SMC. This is again due to the optimal discrepancy for some u being 
greater than e; however, the samples that remain are independent samples. For M/G/l, the results are similar, 
but the ESS/n is lower than the number of discrepancies satisfying e 1, 9% v 12%, indicating that the volume 
of the Jacobians is having a large effect on the variance of the weights. Future work will explore this further. 


noisy. In the M/G/l queuing model, a single server processes arriving customers, which are then 
served within a random time. Customer m arrives at time w m ~ Exp (Os) after customer m — 1, and 
is served in s m (0\,0f) service time. Both w m and s m are unobserved; only the inter-departure 
times Xm are observed. Following mi, we write the simulation algorithm in terms of arrival times 
v m . To simplify the updates, we keep track of the departure times d m . Initially, do = 0 and Vo = 0, 
followed by updates for m > 1: 

Vm = Vjyi— 1 VJ rn Xm = Srri max(0, V m d m — l) dm = d m — 1 H - Vm (22) 

After trying several optimization procedures, we found the most reliable optimizer was simply a 
random walk. The random sources in the problem used for W m (there are M) and for U m (there are 
M), therefore u is dimension 2 M. Typical statistics for this problem are quantiles of x and/or the 
minimum and maximum values; in other words, the vector x is sorted then evenly spaced values for 
the statistics functions / (we used 3 quantiles). The Jacobian is an M x 3 matrix. In our experiments 
<9* = [1.0, 5.0,0.2] 

5 Conclusion 

We have presented Optimization Monte Carlo, a likelihood-free algorithm that, by controlling the 
simulator randomness, transforms traditional ABC inference into a set of optimization procedures. 
By using OMC, scientists can focus attention on finding a useful optimization procedure for their 
simulator, and then use OMC in parallel to generate samples independently. We have shown that 
OMC can also be very efficient, though this will depend on the quality of the optimization pro¬ 
cedure applied to each problem. In our experiments, the simulators were cheap to run, allowing 
Jacobian computations using finite differences. We note that for high-dimensional input spaces and 
expensive simulators, this may be infeasible, solutions include randomized gradient estimates 1(221 
or automatic differentiation (AD) libraries (e.g. Ml Future work will include incorporating AD, 
improving efficiency using Sobol numbers (when the size u is known), incorporating Bayesian opti¬ 
mization, adding partial communication between processes, and inference for expensive simulators 
using gradient-based optimization. 
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